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Abstract. Dark matter (DM) direct detection experiments which are directionally-sensitive 
may be the only method of probing the full velocity distribution function (VDF) of the Galac¬ 
tic DM halo. We present an angular basis for the DM VDF which can be used to parametrise 
the distribution in order to mitigate astrophysical uncertainties in future directional exper¬ 
iments and extract information about the DM halo. This basis consists of discretising the 
VDF in a series of angular bins, with the VDF being only a function of the DM speed v within 
each bin. In contrast to other methods, such as spherical harmonic expansions, the use of 
this basis allows us to guarantee that the resulting VDF is everywhere positive and therefore 
physical. We present a recipe for calculating the event rates corresponding to the discrete 
VDF for an arbitrary number of angular bins N and investigate the discretisation error which 
is introduced in this way. For smooth, Standard Halo Model-like distribution functions, only 
N = 3 angular bins are required to achieve an accuracy of around 10 — 30% in the number 
of events in each bin. Shortly after confirmation of the DM origin of the signal with around 
50 events, this accuracy should be sufficient to allow the discretised velocity distribution to 
be employed reliably. For more extreme VDFs (such as streams), the discretisation error is 
typically much larger, but can be improved with increasing N. This method paves the way 
towards an astrophysics-independent analysis framework for the directional detection of dark 
matter. 
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1 Introduction 

The measurement of the directionality of nuclear recoils in dedicated low-background detec¬ 
tors has long been proposed as a method to detect particle dark matter (DM) [1], The motion 
of the Sun through the Galactic DM halo generates a so-called ‘wind’ of DM particles, which 
would appear to originate from the direction of the constellation of Cygnus. The resulting 
recoil spectrum will be peaked in the opposing direction (or for high mass DM, in a ring 
around this direction [2]). Such a directional signal would be a smoking gun for DM against 
the typically isotropic distribution of residual backgrounds [3-7]. Directional detection may 
also provide a way to circumvent the neutrino floor [8] and to probe the local structure of 
the DM velocity distribution [9, 10]. 

A number of directional experiments are currently in the prototype stage. The pre¬ 
dominant technology is the low-pressure time-projection chamber (TPC), in which nuclear 
recoils leave an 0(1 mm) track of ionised electrons, which can be imaged by drifting the 
electrons (or electron-transport gas) to an anode grid. This technology is currently utilised 
by the DRIFT [11, 12], MIMAC [13, 14], DMTPC [15, 16], NEWAGE [17, 18] and D3 [19] 
collaborations and the problem of reconstructing the recoil direction using this technology has 
been much studied [20, 21]. More recent proposals for directional detectors include nuclear 
emulsions [22], electronic scattering in crystals [23], DNA-based methods [24] and columnar 
recombination in Xenon targets [25, 26]. 

Directional experiments such as these may be the only possibility for probing the full 
3-dimensional velocity distribution of DM, /(v), which is a priori unknown. The standard 
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analysis of direct detection experiments assumes the so-called Standard Halo Model (SHM), a 
simplified model for the DM halo as an isotropic, isothermal sphere of particles [27]. However, 
such a simple description is unlikely to be an accurate description of the Milky Way halo [28- 
32]. In addition, IV-body simulations raise the possibility of more complicated distributions, 
including debris flows [33, 34], tidal streams [35, 36] or a dark disk [37-39]. The impact of 
these astrophysical uncertainties on directional signals has previously been studied [2, 9, 10]. 
While non-standard astrophysics should not severely affect the ‘smoking gun’ directional 
signature of DM recoils [9] , it may still pose a potential problem for future data. In particular, 
trying to reconstruct the DM parameters (such as mass m x and interaction cross section a) 
from a potential signal requires some assumptions to be made about the form of /(v). As 
with non-directional detection techniques, poor assumptions about the astrophysics of DM 
can lead to biased reconstructions of these parameters [40-44] . 

Astrophysics-independent approaches to non-directional experiments have received sig¬ 
nificant attention in the past few years, including so-called ‘halo-independent’ methods [45- 
55] and attempts at parametrising the 1-dimensional speed distribution [40, 42, 43, 56]. There 
have also been several attempts to construct a suitable parametrisation for the 3-dimensional 
velocity distribution, relevant for directional experiments [57-60]. This is significantly more 
difficult than in the non-directional case, owing to the presence of the angular components 
of /(v). A single 1-dimensional function is therefore no longer sufficient to describe the full 
distribution function and a suitable angular basis must be found. Several bases have been 
suggested (e.g. Refs. [59, 60]) but all typically suffer from the same problem, which has not 
previously been addressed in the literature: they allow /(v) to take negative - and therefore 
unphysical - values, which can lead to spurious results. We discuss this problem in more 
detail in Sec. 3. 

In this work, we present a framework which allows the velocity distribution to be 
parametrised in a model-independent way, while guaranteeing that it is everywhere posi¬ 
tive. This consists of an angular discretisation of the velocity distribution, decomposing it 
into a series of 1-dinrensional functions, each of which is constant over a given bin in the 
angular coordinates. This is motivated in the first instance by its ability to describe the sim¬ 
plest directional signal, a forward-backward asymmetry in the scattering rate. However, we 
extend the discretisation to an arbitrary number of bins N in the angular variables, such that 
in the limit of large N, the full 3-dimensional distribution can be recovered. For arbitrary N, 
we also demonstrate how to calculate the Radon transform (the function which encodes the 
angular dependence of the scattering rate) and therefore how to compute the full scattering 
rate. 

Here, we focus on the angular discretisation itself, and do not attempt a full reconstruc¬ 
tion of the DM velocity distribution and particle physics parameters (as in e.g. [40, 42]). 
Instead, we investigate the magnitude of the error introduced by this angular discretisation 
for two different benchmark velocity distributions. To do this, we compare the event rate 
within each angular bin obtained using the full and discrete velocity distributions. Quali¬ 
tatively, we discuss and compare the energy dependence of the event rate within each bin. 
Quantitatively, we can calculate the number of signal events in each angular bin for the full 
and discrete distributions, and compare the result. This allows us to evaluate whether (and 
in which scenarios) fitting to data using such an angular discretisation gives accurate results. 

In order to perform this analysis, we assume that the velocity distribution is known 
exactly. In the analysis of real experiments, a parametrisation of /(v) within each bin must be 
chosen and fit to the data. We leave the choice of a suitable parametrisation of the remaining 
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1 -dimensional functions of v to future work, meaning that the results presented here represent 
a lower limit on the error induced using the discretised velocity distribution. However, we lay 
out a framework which can be used in the analysis of future data to extract coarse-grained 
information about the DM velocity distribution in a way which was not previously possible. 

In Sec. 2, we present the directional detection formalism, including calculation of the 
scattering cross section and the Radon transform, and a brief discussion of the associated 
astrophysical uncertainties. This is followed in Sec. 3 by a discussion of previous attempts 
to overcome these uncertainties. In Sec. 4, we present the discretised velocity distribution 
which is the focus of this work. In Sec. 5, we investigate how this discretised approximation 
compares to the exact result, both at the level of the Radon transform and the event rate. 
Finally, we discuss the significance of these results and plans for future work in Sec. 6 . The 
procedure for calculating the Radon transform from the discretised distribution is presented 
in Appendix B, while we relegate the full derivation of this formula to Appendix A. 


2 Directional formalism 


The directional rate per unit detector mass for recoils of energy Er can be written as [61] 


d R 

dE R dfl q 


±^/c n f\e r )K^A), 


( 2 . 1 ) 


where q = (sin 6 cos 4>, sin 0 sin 0 , cos (/)) is a unit vector pointing in the direction of the recoil. 
We have written the rate R in terms of the DM-proton cross section at zero-momentum 
transfer cr p , which may be spin-independent or spin-dependent; the form factor F 2 (Er), 
which describes the loss of coherence due to the finite size of the nucleus; and an enhancement 
factor Cn, which describes the enhancement of the rate for a nucleus N relative to the DM- 
proton rate. The local DM density is written as po, the DM mass as m x and the reduced 
DM-proton mass as p xp = m x m p / ( m x + m p ). 

For the gaseous targets used in low-pressure TPCs, relatively light nuclei are typically 
used. For example, the DRIFT experiment [11, 12] uses CF 4 as the target gas, with spin- 
dependent (SD) scattering from 19 F nuclei expected to be the dominant interaction with DM 
particles. For SD interactions, the enhancement factor can be written as [62, 63] 


C 8 N D = \(Sp) + a n /a p (S n )\ 2 , ( 2 . 2 ) 

where J is the total nuclear spin, (S p , n ) is the expectation value of the proton and neutron 
spin in the nuclear ground state with maximal magnetic quantum number, and a n /a p is the 
ratio of the DM-neutron and DM-proton couplings. The ratio a n /a p depends on the specific 
model of DM under consideration, although simplifying assumptions are often used, including 
pure-nucleon couplings (a p = 0 or a n = 0) [64] or values motivated by specific models (e.g. 
a n /a p = ±1) [65]. 

The SD form factor can be written in terms of the spin structure functions of the nucleus 
Sij(q), which depend on the momentum transfer q = ^ 2 t 71 nEr, for a nucleus of mass mjv: 


F 2 D (q) = S(q)/S(0), 


with 

S(q) = alS 00 (q ) + a 0 aiS 0 i(<?) + alS u (q) ■ 


(2.3) 

(2.4) 
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The isoscalar and isovector couplings are related to the proton and neutron couplings by 
a -0 = o, p + a n and ai = a p — a n [63]. The spin structure functions (and the proton and 
neutron spin matrix elements) can be calculated using nuclear shell models (see e.g. Refs. [66- 
71]). Here, we focus on Fluorine targets, assuming a p = a n , with spin structure functions 
and matrix elements taken from Ref. [72], The spin structure functions lead to a roughly 
exponential suppression of the recoil spectrum as a function of recoil energy. 

The DM velocity distribution enters into the scattering rate through the Radon trans¬ 
form 

/(umin, q) = [ /(v)d(v • q- u min )d 3 v, (2.5) 

J R 3 

where u m j n is the minimum DM speed required to excite a nuclear recoil of energy Er: 


Dnin — ^min(T'ij) 


m. N E R 


v ^ 


2 

X N 


( 2 . 6 ) 


Geometrically, the Radon transform corresponds to an integral over /(v) on a plane per¬ 
pendicular to q, which has a perpendicular distance u m i n from the origin. Physically, this 
corresponds to integrating over all DM velocities which satisfy the kinematic constraints for 
exciting a nuclear recoil of energy Er in direction q. 

The typical assumption for the form of the DM velocity distribution is the Standard 
Halo Model (SHM). For a spherically symmetric, isothermal DM halo, with density profile 
p(r) oc r~ 2 , the resulting velocity distribution has a Maxwell-Boltzmann form in the reference- 
frame of the Galaxy, [61] 


/(v) 


(27rcj2)3/2 exp 



(2.7) 


with velocity dispersion a v . The corresponding Radon transform also takes the form of a 
Gaussian, 


f(Vi nin,q) = 


(27TO-2) 1 / 2 


exp 


2 


( 2 . 8 ) 


Using the properties of the Radon transform we can obtain the corresponding result in the 
Earth frame by performing a Galilean transformation v —> v + v e , with v e the velocity of 
the Earth with respect to the rest frame of the halo, in which case 


/(^min,q) /(^min + V e • q,q) . (2.9) 

The SHM velocity distribution in the laboratory frame is then 

$/ ~ \ 1 ( (Unin T v e • q) \ 

q) =(^25172-p (—^—)■ < 2 - 10 > 

Within the SHM, one typically assumes values of v e ~ 220 km s -1 and a v = v e /y/2 ~ 
156 km s^ 1 . However, there are still uncertainties on these values of at least 10% [73-75]. 
Introducing a cut off in the distribution at the Galactic escape speed v esc introduces further 
uncertainty into the model [76]. 

Despite its widespread use, the SHM is unlikely to be an accurate representation of the 
DM halo. Observations and N-body simulations indicate that the halo should deviate from a 
1/r 2 profile and may not be spherically symmetric. As a result alternative models have been 
proposed. Speed distributions associated with triaxial halos [28] or with more realistic density 
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profiles [29] have been suggested, as well as analytic parametrisations which should provide 
more realistic behaviour at low and high speeds [30]. Self-consistent distribution functions 
reconstructed from the potential of the Milky Way have also been obtained [31, 32]. 

It is also possible to extract the speed distribution from N-body simulations. Such 
distribution functions tend to peak at lower speeds than the SHM and have a more populated 
high speed tail [77-79]. There are also indications that DM substructure may be significant, 
causing ‘bumps’ in the speed distribution, or that DM which has not completely phase-mixed 
- so-called ‘debris flows’ - may have a contribution [34]. 

Another result obtained from simulations is the possibility of a dark disk. When baryons 
are included in simulations of galaxy formation, this can result in DM subhalos being pref¬ 
erentially dragged into the disk plane where they are tidally stripped [37, 38]. The result is 
a dark disk which corotates with the stellar disk, though with a smaller velocity dispersion, 
typically a® D ~ 50 km s _1 . Such a dark disk could comprise a large fraction (up to 50%) of 
the total DM density, though more recent results [39] suggest a smaller dark disk contribution 
(around 10%), depending on the merger history of the Milky Way. 

Finally, direct detection experiments probe the DM halo on sub-milliparsec scales and 
there is therefore the possibility that DM sub-structures could dominate the local distribution. 
Analyses of N-body simulations suggest that no individual subhalos should dominate the local 
density [80, 81]. However, local structures could still be significant, especially streams of DM 
from the tidal disruption of Milky Way satellites, such as Sagittarius [35, 36, 82], 

Such a wide range of possibilities for the form of /(v) leads to substantial uncertainties 
on the predicted event rate for directional detectors. This in turn can lead to potential bias in 
the reconstruction of other DM parameters, if these uncertainties are not properly accounted 
for. 

3 Problems with parametrising the velocity distribution 

With promising developments in directional detector technology, it is interesting to ask what 
information about the velocity distribution we could, in principle, extract from a direc¬ 
tional signal. Early attempts to address this question involved extending the SHM to allow 
anisotropic velocity dispersions and fitting these new parameters to mock data [57]. The 
possibility of including additional structures, such as streams and dark disks, has also been 
considered [58]. However, such methods are likely to be accurate only if the underlying 
velocity distribution is well-described by the chosen model. 

Alves, Hendri and Wacker [59] investigated the more general possibility of describing 
/(v) in terms of a series of special functions of integrals of motion (energy and angular 
momentum). These can then be fit to data, with around 1000 events required to distinguish 
between the SHM and a Via Lactea II distribution function [83]. However, the special, 
separable form of the velocity distribution requires that the dark matter halo is in equilibrium. 
More troubling is that the resulting velocity distribution is not guaranteed to be everywhere 
positive and therefore not all combinations of parameters correspond to physical distribution 
functions. 

A more general parametrisation for the velocity distribution was recently proposed by 
Lee [60]. In this approach, the velocity distribution is decomposed into products of Fourier- 
Bessel functions and spherical harmonics. This is completely general and does not require 
assumptions about the halo being in equilibrium. Lee also gives an analytic expression 
for the Radon transform of the Fourier-Bessel basis, making this approach computationally 
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efficient. However, this approach can also produce negative-valued, non-physical distribution 
functions, as in the case of the method of Alves et ah. 

In fact, any decomposition in terms of spherical harmonics leads to this problem, because 
the spherical harmonic basis functions can have negative values. It is unclear how this issue 
will affect parameter reconstruction. It may lead to parameter estimates (e.g. for m x or a p ) 
which are unphysical, as they require an unphysical distribution function to fit the data well. 
Without some criteria which determines which coefficients of the spherical harmonics lead 
to strictly positive distribution functions, it may not be possible to reject such parameter 
points. We may attempt to numerically test each parametrised distribution function for 
negative values but for a real function of three parameters /(v) = f(v x ,v y ,v z ) this would 
require a very large number of evaluations, which may not be computationally feasible (and 
does not absolutely guarantee positive-definiteness). In addition, physical distributions may 
occupy only a small fraction of the total space of parameters making parameter sampling 
and reconstruction difficult. 

Even if positive-definiteness could be ensured, it is not clear how to interpret the pa¬ 
rameters reconstructed in this way. Spherical harmonic approximations of typical velocity 
distributions such as the SHM (obtained by integrating out the coefficients of the basis func¬ 
tions) tend to produce distributions which do contain negative values. This is especially true 
when only a small number of basis functions is used. However, the ‘true’ spherical harmonic 
coefficients obtained in this way cannot be obtained by fitting to data (because we would 
reject those which lead to negative values). This may again lead to a bias in the reconstructed 
DM parameters, because the ‘true’ values do not lie in the allowed parameter space. Such 
potential problems have not previously been noted in the literature. 

In order to fit to data, then, it is necessary to decompose /(v) into a series of angular 
components A ®: 

/(v) = f(v,cos6',(j)') = / 1 (u)A 1 (0 , ,^> / ) + f 2 (y)A 2 (6\ ft) + / 3 (u)A 3 (6>', (jj) + .... (3.1) 

We then truncate the series at some order, leaving only a finite number of 1-dimensional 
functions f l (v ) which are unknown. This reduces the problem of attempting to fit a function 
of the 3-dimensional variable v to the problem of parametrising a series of 1-dimensional 
functions, which is much more tractable. Of course, we should be careful that this truncation 
preserves enough angular information to still provide a good approximation to /(v). However, 
as more data becomes available, we can add more terms to the series to capture more angular 
features in the distribution. 

As we have discussed, the spherical harmonic basis may not be an appropriate choice 
for this decomposition. In the next section, I will present an alternative decomposition which 
can guarantee that the velocity distribution is everywhere positive and therefore represents 
a promising and general method for extracting information from directional experiments. 

4 A discretised velocity distribution 

In order to ensure that the velocity distribution is everywhere positive, we propose that the 
velocity distribution be discretised into N angular components: 
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(4.1) 


/(v) = f(v, cos 8', ft) 


f{v) for 0' E [0, n/N] , 
f 2 (v) for O’ E [7r/iV, 27 t/A^] , 

/ fc (u) for O’ E [(fc — 1)7 t/1 V, Aur/iV] , 

f N (v) for 6' E [(IV — l)7r/iV, 7r] . 


Over each bin in /(v) has no angular dependence and depends only on a single function 
of the DM speed. The resulting velocity distribution will be everywhere positive, as long as 
a suitable parametrisation for the f k (v) is chosen which is itself everywhere positive. 

We consider for simplicity only a discretisation in cos O'. We note that in this work we 
will only be considering the azimuthally-averaged event rate (i.e. integrated over the recoil 
angle ft). In this case, we emphasise that no assumptions about the (^'-dependence of /(v) 
are required, as the integral over ft can be performed exactly, and the azimuthally-averaged 
rate depends only on the azimuthally-averaged velocity distribution (see Appendix A). We 
can therefore take the velocity distribution to be independent of ft without loss of generality. 
However, this analysis could be extended to consider bins in the 4> angle, with an additional 
discretisation in ft if required. 

The motivation for this discretised description is that the simplest signal (beyond an 
isotropic N = 1 signal) which can be observed with a directional detector is an asymmetry 
between the event rates in, say, the forward and backward directions. Shortly after the 
confirmation of a dark matter signal at a directional detector, the number of events may still 
be quite small (for example, the roughly 10 events required to distinguish from an isotropic 
background [5], or 30 events required to confirm the peak recoil direction [7]). In this small 
statistics scenario, constraining a large number of free functions is not feasible. However, 
if we discretise /(v) into N = 2 angular components, it may be possible to extract some 
meaningful directional information with only a small number of events. With larger numbers 
of events, N can be increased to allow more directional information to be extracted. 

We show in Fig. 1 some examples of this discretised velocity distribution. We show the 
SHM velocity distribution (top left), as well as the N = 2 (middle left) and N = 3 (bottom 
left) discretised approximations, in all cases integrated over ft. These approximations are 
obtained by averaging the full velocity distribution over each bin in O': 


/T) = 


COS ((k - l)ir/N) - cos(far/A0 J c „ ik ,/ N ) 


f 

J C( 


cos((k— l)7r/N) 


/(v) dcos O'. 


(4.2) 


For comparison, the results for a stream distribution are also shown in the right column. We 
describe these two distribution functions in more detail at the start of Sec. 5. 

Upon discretisation, the distribution which is initially focused in one direction (towards 
O' = 0°) now becomes constant over O' in each bin. For the SHM (left column), there is a 
predominantly forward-going component (cos O' > 0), as well as a smaller backwards compo¬ 
nent (cos O' < 0). For the stream (right column), which has a much narrower dispersion and 
higher peak speed, the N = 2 discretised distribution is almost entirely in the forward direc¬ 
tion. However, we note that due to the discretisation, there is now a sizeable population of 
particles with transverse velocities {O' = 90°), which was not the case in the full distribution. 
In the N = 3 discretisation (bottom row), the velocity distributions become more focused 
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Figure 1. The full velocity distribution (top) as well as N = 2 (middle) and N = 3 (bottom) 
discretised approximations for two examples: the SHM with v e = 220 km s^ 1 , a v = 156 km s -1 (left 
column) and a stream with v e = 500 km s -1 , cr v =20 km s _1 (right column). In each case, we have 
integrated over the cj>' direction and only show f(v,cos9'). The vector v e is aligned along O' = 0. The 
same colour scale is used in each plot. In the case of the full stream distribution (top right), large 
values have been truncated for easier comparison with the other plots. 
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and begin to recover more of the directionality of the full distributions. In the limit of large 
N, the full distribution can be recovered. 

In order to determine the corresponding event spectrum, we must calculate the Radon 
transform from this discretised /(v). We introduce the integrated Radon transform (IRT), 
integrated over the same angular bins as for /(v): 


rcos((j— 1 )tt/N) 

f j (v min) = / / /(nmi^qjdcosfldf/), (4.3) 

J(j )=0 Jcos(jn/N) 

where q = (sin 6 cos 0, sin 6 sin 0, cos q i). 1 In principle, we can define the IRT over a set of bins 
different from those used to define the discretised distribution function. However, using the 
same bins in both cases typically simplifies the calculation of the IRT. A notable exception 
would be to consider fewer bins for the Radon transform than for the distribution function 
(in order reduce possible discretisation error). In this scenario, calculation of the IRT would 
be even simpler. 

The IRT arises from the calculation of the directional recoil rate, integrated over a 
given angular bin, which could then be compared with the number of events observed in 
that bin. We consider the IRT for two reasons. First, it allows us to perform all of the 
relevant angular integrals in the calculation of the Radon transform, eliminating the need 
for computationally-intensive numerical integrals over the angular variables. Secondly, the 
loss of angular information in discretising the velocity distribution means that the full Radon 
transform of this discretised distribution is unlikely to provide a good fit to the data on 
an event-by-event basis. Instead, considering the IRT (and correspondingly binned data) 
should reduce the error which is introduced by using a discretised approximation to the 
velocity distribution. This in turn allows us to parametrise the n-dependence of each angular 
bin and mitigate uncertainties in the velocity distribution. 

The full details of the derivation of (n m i n ) for arbitrary N are included in Appendix. A. 
We include this derivation as a pedagogical tool in the interests of anyone who wishes to mod¬ 
ify or extend the approach presented here (for example, by considering different definitions for 
the bins in cos 6' or cos0). However, for the reader interested in simply calculating / J (n m i n ), 
given a set of f k (v), the algorithm is given in Appendix B. We note that all of the angu¬ 
lar integrals involved can be performed analytically, meaning that computing the discretised 
Radon transform reduces to performing a series of 1-dinrensional integrals over v (one integral 
for each / J ), which can be performed numerically. A python code which implements the full 
calculation of Appendix B is available from the author. 

We note that in the case N = 1, the IRT simply reduces to an integral over all angles, 
exactly recovering the velocity integral which appears in non-directional experiments. Also, 
in the N = 2 case, the IRT takes the following simple form: 


f (Vmin) — 47T 
f 2 (v min) = 47T 



7r / 1 (' u ) + t an 


7 rf 2 (v) + tan 


l 


l 


(A 2 ) 

(A 2 ) 


[f\v) - f 2 (v)] 


dn, 

(4.4) 

dt;, 

(4.5) 


1 The integration over (j) essentially reduces the experiment to what is referred to as 1-dimensional directional 

detection, in which only information about 9 is available. However, it has been shown that this does not 
significantly reduce the discovery power of the detector when compared with 3-dimensional detection [84]. 
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where f3 = v m \ n /v. We have also checked using Monte Carlo calculations that our calculations 
give the correct forms for the IRT in the cases N = 2 and N = 3. In order to do this, we 
randomly sample particles from the discretised distributions. For each particle, we simulate 
a random scattering event and store the values of u m i n and q. We then bin these events in 
angle and compare with the calculated distribution of f J (v m ; n ). The results are in perfect 
agreement with the results expected from Appendix B. 


5 Comparison with exact results 


We now wish to compare these approximate IRTs with the IRTs obtained from the full (non- 
discretised) velocity distribution. To do this, we select a benchmark velocity distribution 
(such as the SHM) and calculate the f k (y ) of Eq. 4.1 by averaging over cos O' in each bin, as 
in Eq. 4.2 and Fig. 1. We then insert these into the algorithm presented in Appendix B to 
obtain the corresponding IRTs. We refer to these as the approximate IRTs. For comparison, 
we use the full Radon transform of Eq. 2.8 to obtain the exact IRTs by integrating over cos 0. 
We fix the peak of the underlying speed distribution to be aligned in the forward direction 
9' = 0°. However, we discuss the consequences for ‘misaligned’ distributions in Sec. 6. 

We discuss qualitatively the differences in shape between the exact and approximate 
IRTs, as any such differences could potentially lead to biases in the reconstruction of parti¬ 
cle physics parameters and the velocity distribution. However, directional direct detection 
experiments do not directly measure the Radon Transform, but rather the differential recoil 
rate dR/dE^dflg of Eq. 2.1. Analogously to the previous section, we define the directionally- 
integrated recoil spectrum: 


dRi f 27r rcos({j-l)-ir/N) 

^-^R 4=0 ■/co.s(j7r/A r ) (1 /l/gll 


d COS 9 d 4> ■ 


The integrated recoil spectrum is then related to the IRT by 

d /?•? 

(X (E R ))F 2 (E R ), 


dEn 


(5.1) 


(5.2) 


where v m i n (Eji) is defined in Eq. 2.6. The form factor F 2 (Er) leads to a roughly exponential 
suppression of the rate (relative to the IRT) with increasing Er. Because of this simple 
proportionality relationship between the IRT and the rate, we do not discuss the shape of 
the rates in detail. 

Instead, we perform a more quantitative comparison between the full and discretised 
event rates. Assuming a particular idealised detector, we compare the expected number of 
events Nj obtained within each bin using the exact and approximate IRTs presented in the 
previous section: 2 




d R? 
d Er 


dEn. 


(5.3) 


In converting from v m i n to recoil energy Er, it is necessary to fix the values of the DM 
and nuclear masses. We therefore take as an example a Fluorine target and DM mass of 
m x = 50 GeV. Direct detection experiments do not probe down to arbitrarily low energies 


“Because the N = 1 case (full average over all angles) is exact, the discretised velocity distribution will 
always recover the correct total number of signal events. We have checked this explicitly up to N = 10 by 
summing the number of events in each angular bin. 
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and typically an energy threshold is set below which either the signal may be dominated by 
unrejected backgrounds or the direction of the recoil cannot be reconstructed. We consider 
a typical threshold energy of 20 keV [85]. 

We fix the normalisation of the event rate by requiring a total of 50 signal events across 
all bins. As previously noted, the anisotropy of a DM signal can be confirmed with around 10 
signal events [5], while the median recoil direction can be confirmed to match that expected 
due to the Earth’s motion with of order 30 events [7]. It is reasonable then to expect that 
the approach presented here might be employed once around 50 signal events are observed. 
In addition, we include a background of 1 event, distributed isotropically. We discuss the 
impact of a larger number of signal events in Sec. 6. 

The error between the exact and approximate results can then be compared with the 
typical Poissonian uncertainty on the expected number of events in each bin A Nj ■M- 
If the discrepancy between the exact and approximate results is smaller than this Poisson 
error in the number of events, we would expect the two methods to be in agreement within 
statistical uncertainties. This means that the data should be approximately as well-fit by 
the discretised velocity distribution as by the full distribution. This in turn indicates that 
the discretised velocity distribution can be used, along with a suitable parametrisation, to fit 
and extract information about the underlying /(v). 

We consider two example distributions, which have already been illustrated in Fig. 1. 
The first is the canonical Standard Halo Model (SHM), with parameters v e = 220 km s^ 1 
and a v = 156 km s _1 . We use this not only because it is so often studied in the literature, 
but also because it has a relatively smooth, simple structure and is not too strongly peaked. 
If the method presented here is to be useful, it should give accurate results for such a simple 
underlying distribution. For simplicity we do not truncate the SHM at the Galactic escape 
speed. The inclusion of such a cut-off would reduce the high-speed DM population, reducing 
the directionality of the signal and therefore improving the results presented here. In the 
reconstruction of real data, an escape speed cut-off could be included easily in the chosen 
parametrisation for f k (y). 

As a comparison, we also consider a stream distribution, modeled using Eq. 2.10, but 
with parameters v e = 500 km s _1 and a v = 20 km s -1 . This leads to a sharply peaked, 
strongly directional distribution function, which allows us to compare the SHM with a more 
extreme case. We note in particular that these parameters lead to a sharper distribution 
(and therefore a more directional rate) than values typically assumed, for example, for the 
Sagittarius stream [35, 36, 82], This stream should therefore be considered a ‘worst-case’ 
scenario which is difficult to approximate accurately. 

5.1 N=2 discretisation 

Figure 2 shows the comparison between the exact and approximate IRTs for the N = 2 
discretisation. The exact IRT is shown as a solid line, while the approximate IRT obtained 
from the discretisation is shown as a dashed line. We show the results for both the SHM 
(blue) and stream distributions (red). Values of u m ; n in the shaded region lie below the 
assumed 20 keV energy threshold for a Fluorine target and 50 GeV DM particle. In the case 
of the SHM, the shape of the Radon transform is well reconstructed in both the forward (left 
panel) and backward (right panel) directions, even with the crude N = 2 discretisation. For 
the forward IRT (j = 1), the ratio between the exact and approximate results is less than 
20%, while for the backward (,j = 2) case, the error is typically larger (50 — 100%), though 
the absolute value is smaller. 
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Figure 2. Exact (solid) and approximate (dashed) integrated Radon transforms, f 3 , defined in 
Eq. 4.3, for N = 2. Results are shown for the SHM (blue) and stream (red) distribution functions. 
The approximate Radon transforms are obtained by discretising the full velocity distribution into N 
angular bins. The vector v e is aligned along O' = 0. The shaded grey region lies below the energy 
threshold of 20 keV for a DM mass of 50 GeV. 


For the stream distribution, the differences between the exact and approximate results 
are more significant. The results using the discretised velocity distribution underestimate 
the forward IRT at low v m i n , while overestimating the backward IRT for all values of 'O’min ■ 
The reason for this can be seen by examining the velocity distributions illustrated in Fig. 1. 
The full stream distribution (top right) is strongly focused in the forward (O' = 0°) direction, 
meaning that particles must scatter through almost 90° to produce a recoil in the backward 
direction. This is kinematically allowed only for a small fraction of the population. Due to the 
angular averaging, however, the N = 2 velocity distribution (middle right) has a significant 
population of DM particles with velocities at right angles (O' = 90°) to the forward direction. 
Thus, the discretised velocity distribution has a greater chance of producing scatters in 
the backward direction. Overall, the discretised distribution is less focused in the forward 
direction than the full distribution, resulting in a reduced asymmetry between the forward 
and backward scattering rates. 

Figure 3 shows the comparison between the number of events Nj in each of the angular 
bins, defined in Eq. 5.3, for the case of N = 2 discretisation. The exact and approximate 
calculations are shown as solid and hatched bars respectively, while the SHM and stream 
distributions are shown separately in the left and right panels. In addition, the ‘exact’ event 
numbers in each bin are assigned an error bar given by the Poissonian standard deviation 

yjy. 

In the case of the SHM (left), the number of events in the forward direction (j = 1) 
is roughly in agreement (within Poisson uncertainties) when calculated using the exact and 
approximate IRTs. However, in the backward direction (j = 2), the approximate calculation 
overestimates the number of events by a factor of 3. As shown in the left hand panel of 
Fig. 2, in the forward direction, the greatest discrepancy between the IRTs occurs at low 
Dnin- When we consider events only above the energy threshold (i.e. above the shaded grey 
region), this discrepancy is therefore minimal. In the right hand panel of Fig. 2, however, the 
fractional discrepancy grows as a function of u m i n , and is greatest above the energy threshold 
of the experiment. This is due to the enhancement in DM particles travelling in the backward 
direction O' = 180° when using the discretised velocity distribution (see Fig. 1) and leads to 
the large discrepancy in the number of backward-going events. 

This problem is even greater for the stream distribution (right), for which neither the 






Bin number, j Bin number, j 


Figure 3. Comparison of the number of events in each angular bin (Eq. 5.3) for N=2 bins, for a 
Fluorine target with 20 keV threshold and m x = 50 GeV. Event numbers are calculated using the 
exact (solid bars) and approximate (hatched bars) IRTs. The SHM distribution is shown in the left 
panel, while the stream distribution is shown on the right. A total of 50 signal events are assumed (as 
well as 1 isotropically distributed background event). The error bars on the ‘exact’ number of events 
are the Poissonian errors A Nj ~ y/Nj- 


forward nor backward event numbers agree within the statistical errors. In the backwards 
direction, the exact calculation predicts much less than 1 event, while the approximate cal¬ 
culation overestimates this by a factor of ~ 30. For such a highly directional distribution, 
such a small number of bins is clearly not suitable. 

We note briefly that in the case of N = 2 discretisation, the approximate calculation 
of event numbers gives results which are almost indistinguishable for the SHM and stream 
distributions. This is because the discretised forms of the SHM and stream distributions 
differ most substantially at low v m i n . This can be seen in the middle row of Fig. 1; eliminat¬ 
ing those particles which are below the threshold speed of ~ 300 km s^ 1 gives two similar 
distributions. Given that there is poor agreement between the number of expected events 
for the N = 2 discretisation, and that the discretised SHM and stream distributions are 
largely indistinguishable, it is clear that such a coarse angular approximation is not suitable 
for fitting to data. 

5.2 N=3 discretisation 

Figure 4 compares the approximate and exact IRTs for the N = 3 discretisation illustrated on 
the bottom row of Fig. 1. In the case of the SHM, the agreement between the two functions is 
improved compared to the N = 2 case. The slope of the event rate is degenerate with the DM 
mass m x , so any error in the shape of the IRTs could translate to a bias in the reconstructed 
value of rn x . In the forward direction (left panel), the fractional error is reduced to at most 
10% and the shape of the IRT appears to be agree closely in all three directional bins. This 
suggests that any such bias should be minimal. 

However, there is still a significant difference between the shapes of the exact and 
approximate IRTs for the stream distribution. For example, in the forward direction (left 
panel), the exact calculation predicts scattering events occurring only with a narrow range 
of Drum from 250 km s -1 to 500 km s^ 1 . This is because almost all particles are moving in 
the forward direction with speed v = 500 km s -1 . The kinematics of the scattering requires 
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Figure 4. As Fig. 2, for N = 3 angular bins. In the right panel, the exact Radon transform for the 
stream is indistinguishable from zero. 



Figure 5. As Fig. 3, for N = 3 angular bins. 


that v • q = , u m in, meaning that scattering in the forward direction (0 = 0°) is only allowed 
for n m ; n ~ 500 km s _1 . Scattering away from the forward direction, but still included in 
the first angular bin (0 < 60°), is kinematically allowed for n m j n > cos(60°) X 500 km s -1 ~ 
250 km s^ 1 . For the discretised velocity distribution, there is a population of DM particles 
initially directed away from the forward direction (0^0°) which can scatter through larger 
angles and still produce scattering events within the first angular bin. These large-angle 
scattering events mean that values of u m i n down to zero are now kinematically allowed. 

Considering the j = 3 angular bin (right panel), we note that the exact stream distri¬ 
bution predicts no scattering events in the backward direction. This is because the particles 
would have to scatter through an angle much larger than 90° to produce events in the j = 3 
bin, which is not kinematically allowed. However, for the discretised distribution (bottom 
right panel of Fig. 1), it is still kinematically possible for particles on the edge of the forward 
bin (0' « 60°) to scatter into the backward direction (0 > 120°). Thus, the approximate 
calculation still predicts a small but non-zero IRT for the j = 3 case. 

Figure 5 shows the event numbers expected in each of the N = 3 angular bins calculated 
using the exact and approximate IRTs. For the SHM (left panel), in the forward direction 
(j = 1), the approximate calculation underestimates the rate by around 13%, while in the 
transverse direction (j = 2), the approximate calculation overestimates the rate by around 
30%. Finally, in the backwards direction (j = 3) the rate is overestimated by a factor of 2. 
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However, these discrepancies are roughly the same magnitude as the typical Poisson error in 
each bin for a sample of 50 events. In contrast to the N = 2 case then, the N = 3 case can 
provide a reasonable approximation to the event rate for the SHM. 

For the stream distribution (right panel), there are significant discrepancies, which re¬ 
main larger than the typical statistical error. In particular, there is a substantial leakage of 
events from the forward direction (j = 1 ) into the transverse direction (j = 2 ) in the approx¬ 
imate calculation. This can be seen in the middle panel of Fig. 4, in which the exact IRT is 
zero above the experimental threshold, while the approximate IRT is substantially non-zero. 
This in turn arises due to the ‘smearing’ of the distribution away from the forward direction, 
as previously discussed. In contrast, in the backward direction (j = 3), the approximate 
calculation now agrees with the exact calculation, giving zero events. As shown in the right 
panel of Fig. 5, the discrepancy between the exact and approximate IRTs arises only at low 
Umin (below the energy threshold of the experiment), where scattering into the backwards 
bin is kinematically allowed. 

5.2.1 The folded rate 

The N = 3 case is important for the scenario where head-tail discrimination is not possible 
when reconstructing recoil tracks. Head-tail discrimination of tracks has previously been 
demonstrated [ 86 ], but may not be possible with 100 % accuracy in near-future detectors 
[21], Within the formalism considered here, detectors without sense discrimination cannot 
distinguish between a recoil with direction cos 6 and another with direction — cos 9. They 
therefore probe the so-called ‘folded’ recoil spectrum: 


d R 


d R 


d R 


dl?ijd| cos 9\ dEjidcos9 + dEjid(—cos9) ^ ^ 

For the N = 2 case, the folded spectrum is entirely isotropic, as the forward and 
backward IRTs differ only in the sign of cos#. However, in the N = 3 case, the transverse 
IRT, given by 


r2n /*l/2 

f 7 (Vmm) = f 2 (v min) = / / /(u min , COS 9) d COS 6 d(j) , 

J<j >=0 J- 1/2 


(5.5) 


is invariant under f(v min .cos 9) —> — cos 9). That is, the transverse event rate ‘folds’ 

back onto itself. Thus, even without sense discrimination, directional experiments will still 
be sensitive to this transverse scattering rate. By comparison, if the forward and backward 
directions cannot be distinguished, the remaining two IRTs (the left and right panels in 
Fig. 4) are folded together, to obtain the longitudinal Radon transform 

r2ir r—1/2 p2n rl 

f L (v min) = / / f (^miri) cos 9 , (f>) d COS 9 d(f> T / / /(u min , cos 9) d cos 9 d(j) . (5.6) 

J<f,=oJ-l J<f>=0 J1/2 

Analogously, we can define the longitudinal and transverse event spectra, dR L,T /dEjt, and 
event numbers N L,T , derived from the corresponding Radon transforms. 

Figure 6 shows the event numbers for this folded rate in the longitudinal and transverse 
bins. For the SHM (left panel), there is a slight improvement in the agreement between 
the exact and approximate results compared with Fig. 5. The backward bin, which was 
previously overestimated by a factor of 2, has been removed and merged with the forward 
bin to form the longitudinal bin. Because the approximate result previously underestimated 
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Figure 6. As Fig. 3, for the folded event rate with N = 3 underlying bins, producing one longitudinal 
and one transverse directional bin (see Eq. 5.4 and associated text). 


the true number of events in the forward bin, this combination of the bins leads to marginally 
improved agreement, though the effect is small in this example. 

In the case of the stream (right panel), the number of events in each bin appears almost 
indistinguishable from those of Fig. 5. This is because the leakage of the events into the 
backwards direction when using the approximate distribution is minimal and therefore the 
effect of folding is also minimal. However, we have demonstrated that this method can still 
be employed when sense discrimination of recoils is not possible and, in general, it should 
lead to a reduction in the associated discretisation error. 

5.3 N=5 discretisation 

We now consider a more finely discretised velocity distribution, namely N = 5. The com¬ 
parison of the exact and approximate IRTs is shown in Fig. 7. The results for the SHM, 
showing typical discrepancies at the 5% level, indicate that the approximate Radon trans¬ 
form does in fact tend to the true Radon transform in the limit of large N. The results 
for the stream again show much poorer agreement. We expect that the agreement will not 
improve significantly until the angular size of each bin is close to the angular extent of the 
underlying distribution function. For the stream, we can see by eye in Fig. 1 that most of 
the distribution is distributed over an angle 9' < 10°, meaning that roughly IV = 18 bins are 
required to prevent the smearing of the distribution function visible for N = 2 and N = 3. 
In spite of this, some structures (such as the peak in the forward direction of the upper left 
panel) are better reconstructed than in the N = 3 case, indicating that additional informa¬ 
tion is still gained by adding more bins. In principle, there is no obstacle to increasing the 
number of bins up to (and beyond) N = 18, as the formalism of Appendix B is applicable 
for arbitrary N. However, the appropriate number of bins would typically be dictated by the 
amount of data available, with a large amount of data required to justify fitting the number 
of parameters associated with N = 18 bins. 

We note that for the bin focused in the backward direction, j = 5 (bottom right), the 
exact and approximate IRTs for the stream are indistinguishable and very close to zero. As 
in the N = 3 case, scattering in the backward direction is not possible for the full stream 
distribution because of the strong directionality. For the discretised distribution, most of 
the distribution is focused in the forward-directed bin (k = 1 in Eq. 4.2). In contrast to 
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Figure 7. Exact (solid) and approximate (dashed) integrated Radon transforms, p defined in 
Eq. 4.3, for N = 5. Results are shown for the SHM (blue) and stream (red) distribution functions. 
The approximate Radon transforms are obtained by discretising the full velocity distribution into N 
angular bins. The vector v e is aligned along O' = 0. In the lower left panel, the exact IRT for the 
stream is indistinguishable from zero, while in the lower right panel, both the approximate and exact 
IRTs are indistinguishable from zero. 


the N = 3 result, however, this is sufficiently directional that scattering into the backward 
direction (into the j = 5 bin) is no longer kinematically allowed, resulting in significantly 
closer agreement in the backward direction. 

In fig. 8, we show the event numbers in each bin for the N = 5 discretisation. For the 
SHM (left panel), there is now close agreement between the exact and approximate results, 
with less than 10% discrepancy for the forward bins j = 1,2 and an error between 20% and 
40% in the transverse and backward bins j = 3,4, 5. For all bins, the agreement between the 
exact and approximate results is significantly smaller than the Poisson uncertainty for the 
example of 50 signal events. This closely matches the expectation from Fig. 7, which shows 
close agreement between the exact and approximate IRTs. 

For the stream (right panel), there is an improved fit between the exact and approximate 
results, with bins j = 1, 2,4, 5 showing agreement within the statistical uncertainty. Leakage 
of events between the two forward bins leads to an slight overestimation of the j = 2 rate. 
However, the biggest problem remains leakage of events into the transverse bin j = 3 which 
gives a significant overestimation of the number of events when using the approximate IRT. 
Once again, this is seen easily in the upper right panel of Fig. 8, where the approximate IRT 
is non-zero above the energy threshold of the experiment. 

We note that while there remains a notable discretisation error for the stream distribu¬ 
tion, it is now possible to distinguish the results for the SHM and stream distribution when 
using the approximate IRTs, in contrast to IV = 2 discretisation. For example, in the j = 3 
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Figure 8. As Fig. 3, for N=5 angular bins. 


bin, the approximate calculation predicts around ~ 8 events for the SHM, compared to ~ 2 
events for the stream. This difference is larger than the error induced by the discretisation 
(a discrepancy of ~ 2 events for both the SHM and stream). This indicates that while the 
N = 5 discretisation may not be able to accurately reproduce the event rate for a stream 
distribution, it may still be useful for distinguishing between different distributions. For ex¬ 
ample, if the N = 5 discretisation were used to fit to data and the results indicated a smaller 
number of events than expected in the j = 3 bin, this could point to deviations from the 
SHM. 

5.4 Forward-backward asymmetry 

Finally, we calculate the forward-backward asymmetry in the number of events, Hfb> and 
compare the exact and approximate results. In terms of the number of forward and backward 
events between the threshold and maximum energies E t h and F? max , 


N f 


N b 



d R 

dE R dtt q 
d R 

dE R dQg 


d cos 9 d 4> dE R 
d cos 9 d(j) d E r , 


the asymmetry is given by 


(5.7) 

(5.8) 


_ N f - N b 
FB N f + N b ' 


(5.9) 


By comparing the values of Afb obtained from the full and discretised distributions, 
we can obtain a measure of how well the discretised distribution captures the directionality 
of the signal, as well as showing how this improves with increasing N. In addition, this 
forward-backward asymmetry would potentially be the first directional signal measured with 
directional detectors and so determining how well it can be recovered is of substantial impor¬ 
tance. Using the current formalism, we can consider only even values of N , for which each 
angular bin lies entirely in either the forward or backward direction. In principle, the calcu¬ 
lation can also be extended to include odd values of N using the framework of Appendix A. 
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Figure 9. Forward-backward asymmetry of the event rate Ape, defined in Eq. 5.9, using the full 
velocity distribution (solid lines) and the discretised velocity distribution (filled squares). Results are 
shown for the SHM and stream distributions in blue and red respectively. Dotted lines show the la 
statistical uncertainties on the measured value of Afb assuming 50 signal events. 


Figure 9 shows the forward-backward asymmetry obtained from the discretised velocity 
distribution for N = 2,4,6,8,10 (filled squares), compared to the true forward-backward 
asymmetry, obtained from the full velocity distribution (solid lines). In the case of the 
stream (red), the event rate is strongly asymmetric (Afb ~ 1 ) as the velocity distribution 
is highly focused in the forward direction. By comparison, the SHM (blue) has a smaller 
asymmetry (Afb ~ 0.9) due to its wider velocity dispersion. Dotted lines show the la- 
statistical uncertainties on the measured value of Afb assuming 50 signal events. 

For both the SHM and stream, the asymmetry is significantly underestimated when 
using only the simple forward-backward (N = 2) discretisation. This matches the analysis 
of the previous sections, which demonstrated that the discretised velocity distribution tends 
to lead to a leakage of events from the forward to the backward direction. However, for both 
distributions, the asymmetry obtained from the discretised distribution rapidly approaches 
the true value with increasing N. For the stream distribution, the approximation converges 
more rapidly than for the SHM. This is because of the strong directionality of the stream 
velocity distribution, so only relatively few bins are required to capture a simple forward- 
backward asymmetry. Even with only 4 angular bins, the error in Afb for both velocity 
distributions is less than 10% and lies within the la statistical uncertainty expected from a 
sample of 50 signal events. 

6 Discussion 

We have demonstrated that for smooth SHM-like distributions, the discretised velocity dis¬ 
tribution allows us to obtain a good approximation to the integrated Radon transform (IRT) 
with relatively few angular bins. Though the N = 2 discretisation is too simple an approxi¬ 
mation, for 3 or more bins the discretised distribution leads to a number of events in each bin 
which matches the true number within statistical uncertainties (assuming 50 signal events). 
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This means that we should be able to parametrise the speed distributions in each angular 
bin f k (v ) and account for astrophysical uncertainties without introducing a large error in the 
number of observed events. For comparison, we have also considered a more extreme and 
highly-directional stream distribution. The disagreement between the true recoil spectrum 
and that obtained from the discretised velocity distribution is significantly larger in this case. 
Though increasing the number of bins N improves the agreement, a rather large number of 
bins (N ~ 18) may be needed to accurately describe the stream distribution. 

The results we have presented here have assumed an idealised future directional detec¬ 
tor. However, application of these results is not necessarily as straightforward in a realistic 
experiment. Conversion from u m ; n to Er requires knowledge of the DM mass rn x and at¬ 
tempting to fit both the velocity distribution and mass simultaneously can lead to spurious 
results [40, 42]. Furthermore, for comparison with experimental data, we must take into ac¬ 
count the fact that experiments have finite angular resolution, typically in the range 20 °-80 °, 
with higher resolution at high energy [21]. Finally, throughout this work, we have assumed 
that the basis for discretising the velocity distribution is aligned with the peak of the under¬ 
lying velocity distribution (i.e. v e aligned along O' = 0). However, this is not known a priori 
and a reliable method of selecting the orientation of the angular basis is required. 

In spite of these open issues, the work presented here is general and conservative. We 
have demonstrated that the stream distribution is more poorly fit than the SHM distribution 
with few angular bins. However, the stream is an extreme example, as discussed in Sec. 5, 
and it is unlikely the discretisation error will be any greater than observed in that case. 
Other possibilities for the velocity distribution would tend to be less directional than the 
stream example, leading to better agreement between the exact and approximate rates. For 
a dark disk, the low value of v e means that the velocity distribution appears almost isotropic 
in the lab frame, which will only reduce the angular discretisation error. A misalignment 
between v e and O' = 0° would also lead to a distribution more isotropic in O '. Similarly, 
including finite angular resolution will smear the velocity distribution (in a style similar 
to the discretisation), reducing the directional asymmetry and therefore the discretisation 
error. We have also demonstrated that this method is applicable in detectors where sense 
recognition is not possible. 

For a realistic analysis, all of these issues should be taken into account. Based on mock 
(or future) data, the optimal direction for O' = 0° should be selected (potentially based on 
the observed median recoil direction [5, 7]). An appropriate number of bins N should then be 
chosen, which may be influenced by the angular resolution of the experiment and the number 
of observed events. In some scenarios, such as for strongly peaked signals, deviations from 
the smooth SHM may be obvious from the data. In this case, a definition of the angular bins 
different from that used here (including differently-sized bins) may be optimal to reduce the 
discretisation error. Using the results of Appendix A, such a scenario can be straightforwardly 
accommodated, although we leave the issue of choosing the number and size of the bins to 
future studies. 

In order to make more quantitative statements, we have focused on a specific scenario, 
in which 50 signal events (and 1 background event) are observed. This has allowed us to 
compare the typical Poisson uncertainties on the number of observed events with the error 
induced by using the discretised velocity distribution. Though a signal consisting of 50 events 
is somewhat arbitrary, it is representative of the lower limit on the number of events required 
before such an analysis would be reasonable. As the number of signal events increases, the 
typical statistical uncertainty decreases. Eventually, this uncertainty becomes smaller than 
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the error between the approximate and exact results. Increasing the number of angular bins 
would then reduce the discretisation error and increase the statistical uncertainty on the 
number of events in each bin, reducing any bias which may be induced by the discretisation. 
This therefore provides a natural scheme for deciding how many angular bins are required, 
based on the number of signal events observed. 

Once the number of bins N has been fixed and a suitable parametrisation for the radial 
functions f k {y) has been chosen, the full parameter space (of both DM particle physics 
parameters and astrophysics parameters) should be fit to the data. We note in particular 
that by summing the radial functions f k (v), the 1-dimensional speed distribution is recovered, 
meaning that non-directional experiments can easily be included in this framework. As 
emphasised in this section, the fitting process will be highly dependent on the underlying 
DM and experimental parameters, as well as on the chosen parametrisation for f k (v). We 
have therefore focused here on the discretisation of /(v) and left more involved investigations 
for future work. 


7 Conclusions 


In this work, we have presented an angular basis which can be used to parametrise the DM 
velocity distribution. This involves discretising the velocity distribution into N angular bins: 


/(v) = f(v,cosd',<f/) 


/ 1 (u) for O' € [0, n/N] , 
f 2 (v) for O' e [tt/N,2tt/N] , 

f k (v) for O' € [(k — 1 )tt/N, kn/N] , 

f N (v) for 0' G [(IV — l)7r/IV, 7r] . 


(7.1) 


In Appendices A and B, we have provided recipes for calculating the corresponding directionally- 
integrated Radon transforms (IRTs), which appear in the directional event rate. Alternative 
methods for parametrising the DM velocity distribution (such as Refs. [59, 60], and in par¬ 
ticular those based on spherical harmonics) may lead to negative values /(v), leading to 
unphysical distribution functions. It is not clear what affect this might have on parameter 
reconstruction or how this problem can be mitigated. The advantage of the basis presented 
here is that it guarantees that the resulting distribution function is everywhere positive and 
therefore physical. 

We have investigated the possible size of discretisation errors in the directional event 
rate by comparing the IRT in each bin obtained from the discretised and full velocity dis¬ 
tributions. The simplest possible discretisation (N = 2) is unsuitable for fitting to data, as 
it substantially underestimates the forward scattering rate, while underestimating the back¬ 
wards scattering rate. However, starting from N = 3, the standardly-assumed SHM can be 
well fit by the discretised velocity distribution. We have shown that once the DM nature 
of the signal is confirmed (with around 50 events), the N = 3 discretisation allows us to 
calculate the number of events in each angular bin and obtain good agreement with the true 
event numbers (within statistical uncertainties). Increasing the number of bins to IV = 5 
improves this agreement, giving discretisation errors in the range 10-50% depending which 
angular bin is considered. 


- 21 



For comparison, we also consider an extreme, highly directional stream distribution. In 
this case, the discretisation error is much larger and exceeds the statistical uncertainties in 
the N = 5 case. An estimated N = 18 angular bins would be required to accurately map 
the full velocity distribution. However, increasing N does reduce the discrepancy between 


the exact and approximate results. In addition, with as few as N = 5 bins, it is possible 


to distinguish between the SHM and stream distribution. This suggests that even though 
the discretised distribution gives a poor approximation to the stream, it may be useful in 
distinguishing different distributions, even for relatively few bins. 

Finally, we have considered the forward-backward asymmetry Afb hi the event rate, 
calculated using both the full velocity distribution and the discretised velocity distribution. 
The latter calculation rapidly converges to the correct value with increasing N. For both the 
SHM and stream distributions, the error in Ape is below 10% for N > 4 bins, within the 
statistical error expected for 50 signal events. We have also demonstrated that the method 
presented here can still be used when head-tail discrimination of the recoil tracks is not 
possible, in which case N > 3 bins are required. 

In this work, we have only considered discretising the angular component of the DM 
velocity distribution, leaving the N radial functions f k (v) fixed to their ‘true’ values. In future 
work, it will be necessary to combine the discretisation presented here with a parametrisation 
of these radial functions (such as the parametrisations of Refs. [40, 42, 44]), in order to 
reconstruct the velocity distribution (and other DM parameters) from mock data. In addition, 
several questions are still to be addressed, including how to decide on the optimal ‘forward’ 
direction for the discretisation, and how realistic angular resolution would impact the results 
presented here. However, these additional considerations are only expected to improve the 
agreement between the true and approximate event rates. This work represents an initial step 
towards a parametrisation of the full DM velocity distribution, which would allow future data 
from directional experiments to be analysed without astrophysical uncertainties and which 
would potentially allow the velocity distribution itself to be probed and reconstructed. 
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A Calculating the Radon transform 

In this appendix, we derive the Radon transform corresponding to the discretised velocity 
distribution. The final result can be found in Appendix B. The Radon transform of /(v) is 


defined by: 



) d 3 


V . 


(A.l) 


We write the coordinates of the velocity and recoil momentum as 
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(A.2) 


v = v (sin O' cos ft, sin O' sin ft, cos O') 
q = (sin 0 cos <f), sin 0 sin (j), cos 0) , 

and consider here only the azimuthally integrated Radon transform. We write this as 
f(v min ,cos 0), which is given by 


/fanimCOsfl) = / f (Wmin, COS 0, (f>) d(f> 

Jo 

= J (^ j 3 /(v) S (v • q - U min ) d 3 v^ d (j) 

= j 3 /(v) 5 (v • q - u min ) d<^ d 3 v 

= [ f(v)D(v m i n ,cosO,v) d 3 v . 

J R 3 

We expand the <5—function explicitly in terms of the angular coordinates: 


(A.3) 


(5 (v • q — U min ) = -5 (sin 0 sin O' cos {d> — ft) + cos# cos# 7 — v m \ n /v) = -5 [gift]) 
v v 

We rewrite the argument of the 5— function as a function ft. 

x ( (x\\ ^ ~~ 4>i) 

Here, we sum over those values of ft satisfying g(ft) = 0: 

, /3 — cos 0 cos O' 

cos{ft -(p) = - ■ Q ■ a , — = « , 

sm 0 sm O' 

where we have also defined [3 = v m \ n lv. The solutions for cj> E [0, 27r] are: 


(A.4) 


(A.5) 


(A.6) 


0i = ft + cos -1 a, 

for ft E [0, 27T 

4>2 = ft + 2vr — cos -1 cc, 

for ft E [0, cos 

4*3 = ft + cos -1 a — 27r, 

for ft E [27T — 

04 = (f) — COS Cl, 

for ft E [cos -1 


2tt] . (A.7) 

We note that these solutions exist only for f3 E [0,1] (or equivalently v > w m i n ) and for 
a E [—1,1], otherwise Eq. A.6 cannot be satisfied. If these constraints are satisfied, there 
exist exactly 2 solutions for a given value of ft and therefore 2 5-functions in Eq. A.5. 

For the derivative of g( ft) we obtain 


g\ ft) = — sin 0 sin O' sin(</> — ft). 
Substituting the values of 01,2,3,4-. we see that 

\g' (0i ,2,3,4 ) | = \j (sin 0 sin O') 2 — (f3 — cos 0 cos O') 2 . 


(A.8) 


(A.9) 
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Each of the two 5-functions therefore contributes the same amount to the integral, regardless 
of the value of <p'. We can now perform the integral over (j>: 


1 f Z7T 

D (u min , cos 9, v) = - / 5 d 

V Jo 

2 C(a) 


v\J (sin 9 sin 9') 2 — (j3 — cos 9 cos 9')' 

= -C(a)I(/3, cos 9, cos 9')@(v — u m i n ) 
v 


-.0(v- v min ) 


(A.10) 


where C(a) = 1 for a € [—1,1] and vanishes otherwise. We note that D is independent of 
the azimuthal angle <J>. The constraint a € [—1,1] is satisfied for cos 9' € [x_,a:+], where 


x± = P cos 9 ± \/l — /3 2 sin 6 . 


We therefore obtain 


rx-\- /*oo 

f fu min , cos 9 ) = 2 / / /(u, cos 9')I(/3, cos 0, cos 0 7 ) u du dcos 9 ', 

J X— J Vxnin 

where we have performed the (j)' integral over the velocity distribution, and defined 

r2n 


/*Z7T 

f(v,cos9')= f(v, cos 9\ J)') dcj )', 
JO 


(A.ll) 

(A.12) 

(A.13) 


because I does not depend on J>. We note from this that the azimuthally-integrated Radon 
transform is unaffected by the exact (//-dependence of /(v). Instead, the azimuthally- 
integrated Radon transform depends only on the azimuthally-integrated velocity distribu¬ 
tion. For the framework considered here, then, we can assume that /(v) is independent of 
4>' without loss of generality. 

We will consider a velocity distribution which is discretised into N angular pieces: 


/(v) = f (v, cos 9\ <j/) = / 


f\v) 

for 9' € [0,7t/AT] , 

f 2 (v) 

for 9’ € [n/N,2ir/N\ , 

f k (v) 

for 9’ G [(A; — 1)tt/N, kn/N] 

J N (v) 

for 9’ G [(IV — l)7r/lV, 7r] . 


(A.14) 


We would then ultimately like to calculate the directionally-integrated Radon transform f 3 , 
integrated over the same bins in cos 9 as for the velocity distribution, 


rcos((j-l)n/N) ^ 

f j [v min) = / /min , COS 9) d COS 6 

J cos(jn/N) 


= 2 


f*CL * ^ roo / /*x \ 

/ / ( / /(^) cos 9')I((3, cos 9, cos 9') dcos 6' j udr>dcos$. 

Jan Jv min \j X- J 


(A.15) 
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Here, we have introduced the notation a k = cos(kn/N) and defined f(v, cos O') = 27 xf{v, cos O', (j)'). 
We divide the integration range for cos O' into different regions: 

cos O' E [x-, afc__i] 
cos O' E [a k _-i,a k _- 2 ] 

: (A.16) 

cos 0 E [u/c++i>®fc+] 
cos O' E [a k+ , x+], 

where k± are dehned such that: 

x± E [a k ±, Ufc±-i] • (A. 17) 

These regions are chosen such that f(v,cosO') is independent of cos O' within each region. 

We note that if k+ = k-, there is only one region: cos O' E [x-,x+\. Using these definitions 
we can therefore rewrite the term in brackets in Eq. A. 15 as 


f f{v, cos 0')I(/3, cos 0, cos O') d cos O' 

J X- 


PO'k_ — 1 

= 2-Kf k -(y) / I(fi, cos 0, cos O')& cos O' 

J X- 

fc + + 1 ra . k -1 

+ ^ 27t f k {v) / I(f3, cos 0, cos O') d cos O’ 


k=k— — l 
ck. 


r x + 


+ 2nf k +(v) / I((3, cos 0, cos 9') d cos O'. 


dk_L 


We can now explicitly perform the integral over cos# 7 , 


r a k -1 


I ((3, cos 0, cos O') d cos O' 


" a k 


rak-1 


d cos O' 


'CLk 


■sj (sin 0 sin O') 2 — (/3 — cos 0 cos O') 2 


_i ( cos O' — j3 cos 0 
sm 1 


— (3 2 sin 0 




- Clfc 


= F(a k ~ i, cos 0) — F(a k , cos 0 ), 


where we have defined 


We note that 


F{y, cos 0) = sin 


. _i j y — f3 cos 0 


F(x±, cosO) = ±- , 


\Jl — (3 2 sin 0 ) 


7r 


(A.18) 


(A.19) 


(A.20) 
(A.21) 
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for all cos 6. We also note that F(y, cos 6) can be integrated analytically: 


F(y , cos 0) cos 0 


= x sin | 

+ y taiW 1 
1 

+ - tan 

1 _i 

-tan 


y - fa \ 

Vl — x 2 -\/l — (3 2 y 
x - y/3\ 

. yft ) 

- y 2 - /3 2 - x + y(3( 1 + x) 

v (y - P)Vt 

— y 2 — (3 2 + x — y(3( 1 — x) 

v {y + /3)Vt 


(A.22) 


= J(x,y)+ C 

where it is understood that J(x, y) is a function of /3 and where 

t = (1 - x 2 )(l - /3 2 ) - (y - (3x) 2 . (A.23) 

It is now possible to also complete the integration over cos0. Assuming that k± do not 
change over a range cos 0 6 [bi, bi + 1 ], we can combine Eqs. A.18 and A.22 to obtain 


H+i / rx+ 


f(v, cos 6')I((3, cos 9, cos O') dcos^J dcos@ 

= 2irf k ~ (v) (r^iPi+i ~ bi) + J(b i+1 ,a k _-i) - J(bi, a k _-i)"j 

k ++1 

+ ^2 27T f k ( v ) {J{bi+h a k-i) - J(bi,a k _i) - J(b i+ i,a k ) + J{b h a k )) 

k=k— — l 

+ 2i Tf k +(v) (^{bi +1 ~ bi) - J(b i+1 ,a k+ ) - J(6i,a fc+ )^) . 


(A.24) 


However, we must take into account the fact that the values of x± and therefore of k± depend 
on (3 and cos 9. We now discuss how to divide up the integration regions of v and cos 6 such 
that we can apply Eq. A.24. 

As an illustration, we show in Fig. 10 the values of x + (x_) as solid (dashed) lines as a 
function of cos0, for a fixed value of /3. In evaluating Eq. A. 15, we wish to integrate over the 
region enclosed by the solid and dashed lines, for the relevant range of cos 6. As an example, 
we show with horizontal and vertical dotted lines the edges of the discretised regions in cos 6 
and cos 9' for N = 3. That is, the dotted lines show the values cos(n7r/lV) for n = 0,1, 2, 3. 
If we wish to evaluate f 1 (v m ; n ), we need to integrate over the shaded region in Fig. 10. 

In Fig. 11, we show x± for different values of (3. We give two examples, for two different 
ranges of f3\ f3 < 1/2 (black) and (3 > \/3/2 (blue). The values of k± (and therefore the 
relevant f k (v)) will clearly depend on the value of (3, so we must be careful to account for 
this properly. 

Calculation of the values of k± and the corresponding bin edges in cos 6 ({6j}) is involved, 
though not technically difficult. We therefore sketch the procedure here. We begin with the 
definition of x + : 
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Figure 10. Integration limits for Eq. A.15 as a function of cos (9, for a fixed value of /3. The limits 
x+ and X- are shown as solid and dashed lines respectively. The dotted vertical and horizontal lines 
mark the values cos(nn/N) for n = 0,1,2,3 and N = 3. If we wish to calculate / 1 (t' m in) J we must 
perform the angular integral over the shaded region. 



cos 6 


Figure 11. As Fig. 10, but for two different values of j3. The black lines show x± for a value of /3 in 
the range /3 < 1/2 while blue lines show x± for a value in the range /3 > y/3/2. 


x + = x + (cos 0) = f3 cos 9 + \Jl — f5 2 sin 6 . (A.25) 

By straight-forward differentiation, we observe that this function has a maximum at cos 6 = f3. 
For concreteness, we consider (3 in the range, 

(3 £ [cos(nir/N),cos((n — l)ir/N)\ for n = 1, 2,iV/2 . (A.26) 

We must consider three separate cases, depending on j (i.e. which cos 6 bin we are consider¬ 
ing): 

(i) j <{n- 1) 

In this case, cos 6 > f3, meaning that dx + /dcosO < 0, so .t_|_ is monotonically decreasing 
in cos 9. 
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(ii) j > (n + 1) 

In this case, cos 9 < (3, meaning that dx+/dcos 9 > 0, so is monotonically increasing 
in cos 9. 

(hi) j = n 

In this case, the maximum of x+ lies in the j th bin. 

For each value of j, we can now determine the maximum and minimum value of x+ 
within that cos 9 bin by using the definition in Eq. A.25 . This allows us to determine the 
value of (i.e. the cos 9' bin into which x+ falls). In fact, it can be shown that for each 
value of j (apart from j = n ), x + crosses one of the cos 9' bin edges exactly once, when 
cos 9 = cos 9 J + , satisfying 


x_l_(cos6^) = cos(lir/N ), (A.27) 

for some l = 1,2,..., N. This means that for a given j, there are two distinct regions in cos0 
(cos 9 < cos 9 J + and cos 9 > cos 9 J + ) each with a different value of k + . 

We can perform a similar analysis for x_ to obtain the values of for a given n and j , 
as well as the values of cos 9 ] _, where the curve x_ crosses from one bin to the next. However, 
we can only apply Eq. A.24 if both k + and k- do not change over the range cos 9 £ [bi, bi + 1 ]. 
It is clear then that we must subdivide each cos 9 bin into 3 smaller bins. The edges of these 
sub-bins will be cos0^ and cos$/_, at which the value of either or k- changes. 

Finally, we need to determine which of these crossings cos 9± occurs first. The crossings 
occur at the same value of cos 9 if cos 9+ = cos #/_, which occurs for (3 = cos ((n — l/2)7r/Af). 
For each value of n, we therefore need to distinguish two regimes, depending on which of the 
crossings cos 9± is larger. It is helpful then to further subdivide the range of /3, writing 

(3 £ [cos(m7r/(2iV)), cos((m — l)vr/(27V))] for m = 1,2,..., IV, (A.28) 


from which we obtain 


n = 


(m + l)/2 
m/2 


for m odd 
for m even. 


(A.29) 


We can now determine which of cos 9± is greater, depending on whether m is even or odd. For 
a given value of m then, the values of the sub-bin edges in cos 9 ({&i}) are now well determined, 
along with the values of k± within each sub-bin. This allows us to apply Eq. A.24 to each of 
these sub-bins. 

Finally, we can perform the integral over v described in Eq. A. 15. However, we note 
that we should sub-divide the range of integration, depending on the value of m. That is, 
the integral over v becomes, 




Vmin/ a m/2 
min/^(m — 1)/2 


dx, 


(A.30) 


where we remind the reader that = cos(fc-7r/./V). For each term in this sum, the value of 
m does not change within the range of integration, so we can apply the results which have 
been obtained so far. We give the full expression for the IRT, as well as the values of k± and 
cos (calculated as described here) in Appendix B. 
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B The directionally-integrated Radon transform (IRT) 

The discretised velocity distribution, with N angular bins, is defined as 


/(v) = f(v, cos O', 4>') 


/ 1 (n) for O' € [0,ir/N] , 
f 2 (v) for O' e [7r/iV, 27 t/A^] , 

f k (v ) for 6' € [(A; — l)n/N, kir/N] , 

f N (v) for O' € [(N — l)vr/iV, 7r] . 


(B.l) 


The integrated Radon transform (IRT), integrated over the same angular bins, is defined as 


f j (y min) 



/(^min,q)dC0s6»d^, 


(B.2) 


where have defined the shorthand a\~ = cos(kir/N). The IRT is given in terms of the distri¬ 
bution functions f k (v ) as 


= 47T 

N 

E 

rV min/ a 

m/2 

v dv 

E 

|(/ t +W + /‘ 1 (»))(i> i +i 


m= 1 ' 

' v min/ a (m— 1)/2 

=0.1,2 

- z 

+ (1- 

_ $k\ 

ki)f k ~ 

(v) ( J(bi+1 : 

> a k i _- 

1) ~ 

- fi¬ 

- $k\ 

kOf k + 

(v) (j(bi+1. 

• a k\) 

- J( 6 „a^)) 

fe 1 

-1 






+ ^2 / fe (^) (J(b i+ i,a k -i) - J(bi,a k ~i) 

k=k l + +1 

J(pi+ 1 ; Qfe) T ®fc)) ] • 


The function J(x,y) depends implicitly on the speed v and is defined below. For each value 
of j, we integrate over N ranges in v, labeled by the index m. We also sum over three ‘bins’, 
labeled by the index i = 0 , 1 , 2 , which are distinguished by their values of the integers k±. 
We note that the values of k± will depend on the values of j and m under consideration, as 
well as the index i. We write k-t = (&;±, k±) and give the values below. 

For odd m, we define n = [m + l)/2 and bin edges 


&o = cos (jn/N) 

b\ = cos (7 + (j + n — 1)7 t/N) 

b -2 = cos (7 + (j — n)7r/N) 

63 = cos ((j - 1 ) 7 t/N) , 

with 7 = cos -1 (u min /u). The corresponding k-t values are: 


(B.4) 
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and 


k+ 


' {n — j,n — j, n-j + 1) 

< ( 1 , 1 , 1 ) 

AJ ~n + l,j -n + l,j -n) 


for j < (n — 1 ) 

for j = n 

for j > (n + 1 ), 


(B. 5 ) 


'(n + j,n + j - l,n + j - 1) 

< (N,N,N) 

k (2TV — n — j + 1, 2TV — n — j + 2, 2N — n — j + 2) 


for j < (TV — n) 
for j = (TV + 1 — n) 
for j > (TV + 2 — n). 


(B. 6 ) 


For even m, we write n = m/2 and bin edges 


bo = cos ( jn/N) 

b\ = cos (7 + (j — 77 ) 77 /TV) 

62 = cos (7 + (j + n — l)n/N) 

63 = cos ((j - 1 )7 t/TV) . 

The corresponding k-t values are: 


(B. 7 ) 


k+ 


' (n — j, 77 — j + 1,77 - j + 1) 

< ( 1 , 1 , 1 ) 

Sj - n+ 1’i - 77, j - 77) 


for j < (77 — 1 ) 

for j = n 

for j > (n + 1 ), 


(B. 8 ) 


and 


k_ 


' (n + j,n + j, 77 + j - 1 ) 

< (iV,TV,TV) 

^ (2TV — 77 — j + 1 , 2TV — n — j + 1, 2TV — 77 — j + 2) 


for j < (TV — 77) 
for j = (TV + 1 — 77) 
for j > (TV + 2 — 77). 


(B. 9 ) 


Finally, the angular function J is given by: 


J(x,y) 


= x sin 
+ y tan ^ 1 


y - fix 


y/\ — x 2 \J\ — (5 2 / 

f x — y/3' 


Vt 


, f t -1 fl-y 2 - P 2 -x + y(3(l + x) 
H— tan --- 

2 V (y- P)Vt 

l t -1 /^l - y 2 - + x-7/^(1-x) 

2 V (y + P)Vt 


(B. 10 ) 


where we define /3 = v m \ n /v and t = (1 — x 2 )(l — (3 2 ) — (y — (3x) 2 . We note that the angular 
integrals in the Radon transform lead to this analytic form for arbitrary TV, meaning that 
only 1 -dimensional integrals over v need to be performed. A python implementation of this 
algorithm is available on request from the author. 
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